import java
from Matrix import Matrix
from affymetrix.fusion.psi import *;
from affymetrix.fusion.cel import *;
from affymetrix.fusion.chp import *;
from affymetrix.fusion.cdf import *
from affymetrix.fusion import *;
from java.util import *;
import shelve;

#import org.systemsbiology.inference
#from Matrix import nonNull
# package to support affy data processing.
#### RMA Method ####
# background adjustment
# quantile normalization <-- (from ICSB Seattle)
# summarization (median polish)
# map probe sets to gene identifiers.
# then provide a list, or shorted cel file that includes only the genes represented on both chips.


class AffxCDF(java.lang.Object):
	cdf=None
	def __init__(self,filename=None):
		r = self.loadFrmFile(filename) # see if we can load the serialized cdf file (should be faster)
		if False==r:
			r = self.loadFrmWeb(filename)
			if False==r:
				print "error loading CDF from File or Local Cached Database", filename

	def loadFrmWeb(self,filename):
		return False

	def loadFrmFile(self,filename):
		self.cdf = FusionCDFData();
		self.cdf.setFileName(filename);
		if self.cdf.exists() == False:
			print "The CDF file does not exist."
			return False;
		# can I attach a progress monitor here..
		if self.cdf.read() == False:
			self.cdf = None
			print "Failed to read the CDF file.";
			return False;
		#read probe sets
		print "Num ProbeSets:", self.cdf.getHeader().getNumProbeSets()
		print "Array type = " + self.cdf.getChipType();

	def processCelFiles( F ):
		N = self.cdf.getHeader().getNumProbeSets()
		info=FusionCDFProbeSetInformation()
		grp = FusionCDFProbeGroupInformation()
		cell = FusionCDFProbeInformation()
		for idx in range(0,N):
			self.cdf.getProbeSetInformation(idx, info)
			for g in range(0,info.getNumGroups()):
				info.getGroup(g,grp)
				C = grp.getNumCells()
				for c in range(0,C):
					grp.getCell(c, cell)
					for cell in F
					# F list # now process all probesets in each CEL file
					# 
					x=cell.getX()
					y=cell.getY()
					
					cell.clear()
		### get probe indexes mapped to group.
			grp.clear()
			info.clear();
		return True

def kNN(E):
	X = nonNull(E)
	X = map( lambda x: x(1), X)
	[R,C]=size(E)
	ex = Set(range(1,C))
	ex.setdifference(X)
	for i in R:
		m=nonNull(E[:,ex])
		n=nonNull(E[:,x])
		q=m.intersect(n)
		[u,s,vt]=svd(m[q,:])
		# then use the s*vt[ex,ex]' svd decomposition to project for E[i,ex]

def importCHP(chpfile):
	FusionCHPLegacyData.registerReader();
        chp = FusionCHPDataReg.read(chpfile);
        if (chp == None):
		print "Failed to read the CHP file."
		return;

	fusionchp = FusionCHPLegacyData.fromBase(chp);
        print "Alg name = " + fusionchp.getHeader().getAlgName();
        System.out.println("Alg version = " + fusionchp.getHeader().getAlgVersion());
        System.out.println("Array type = " + fusionchp.getHeader().getChipType());
        System.out.println("# probe sets = " + fusionchp.getHeader().getNumProbeSets());
        params = fusionchp.getHeader().getAlgorithmParameters();
        for i in range(0,params.size()):
		param = params.elementAt(i)
		System.out.println(param.getTag() + " = " + param.getValue());
        params = fusionchp.getHeader().getSummaryParameters();
        for i in range(0,params.size()):
            param = params.elementAt(i);
            System.out.println(pFusionCELDataaram.getTag() + " = " + param.getValue());
        # maps probes to probesets
	
        System.out.println("RawQ from alg params = " + fusionchp.getHeader().getAlgorithmParameter("RawQ"));
        System.out.println("RawQ from summary params = " + fusionchp.getHeader().getSummaryParameter("RawQ"));

def importCEL(filename):
	cel = FusionCELData();
	cel.setFileName(filename);
	if cel.exists() == False:
		print "The CEL file does not exist."
		return;
		
	if cel.read() == False:
	        print "Failed to read the CEL file.";
		return;
	print cel.getHeader()
        print "Algorithm name = " + cel.getAlg();
        print "Array type = " + cel.getChipType();
	n = cel.getCells();
        sum = 0;
        for i in range(0,n):
		sum += cel.getIntensity(i);
		
	# read in CDF file.
	# cell files map to probe sets. PM/MM   # need to write the process handler here.
	# for each probe set compute avereage and take background correction
        print "The average intensity is " + (sum/n).__repr__()


def normalizeByGene(cell,dictionary):
	pass

def getBkgLevel(affy):
	# what algorithm do I use?
	pass

def qNorm(affy):
	pass

def summarize(affy):
	pass

def logratios(affy1,affy2):	
	n1 = getBkgNoiseLevel(affy1) #if the data is RMA then we can skip this step  # can we use the SVD?
	n2 = getBkgNoiseLevel(affy2)
	qn1 = quantNorm(n1)
	qn2 = quantNorm(n2)
	c1 = summarize(qn1)
	c2 = summarize(qn2)
	
	if c1 > c2:
		R = (c1 - c2) / n1	
	else:
		R = (c2 - c1) / n2
	return R

